% irf.m

function z = irf(varargin)
% Syntax
% =======
% 
% function response = irf(PI,PHIs,PHIn,...
%     T2zq,T2z,T2qq,T2qv,T2q,...
%     T3zq,T3z,T3qq,T3qv,T3q,...)
%
% Options
% ========
% 
% -  'state=': initial state vector (default: zeros)
% -  'shock=': initial shock impulse vector (default: zeros). Can be a
% vector with the appropriate number of shocks or an array in which case
% the array is split into vectors of shocks to be used in sequential
% periods. 
% -  'T=': number of periods to simulate (default: 40)
% 
% Description
% ============
% 
% Computes non-linear impulse responses and/or transitional dynamics based
% on Johnson, King and Lie (2014).
%
% period 1: initial state
% period 2: shock hits
% 
% Assumes a state space representation as in the appendices to
% Johnston, King, and Lie (2014). 



[state,varargin] = parseOption(varargin,'state',[]) ;
[shock,varargin] = parseOption(varargin,'shock',[]) ;
[T,varargin] = parseOption(varargin,'T',40) ;

if numel(varargin)>9
    order = 3 ;
elseif numel(varargin)>3
    order = 2 ;
else
    order = 1 ;
end

% first order input arguments
PI = varargin{1} ;
PHIs = varargin{2} ;
PHIn = varargin{3} ;
if order>1
    % second order input arguments
    T2zq = varargin{4} ;
    T2z = varargin{5} ;
    T2qq = varargin{6} ;
    T2qv = varargin{7} ;
    T2q = varargin{8} ;
    nq2 = size(T2zq,2) ;
    if order>2
        % third order input arguments
        T3zq = varargin{9} ;
        T3z = varargin{10} ;
        T3qq = varargin{11} ;
        T3qv = varargin{12} ;
        T3q = varargin{13} ;
        nq3 = size(T3zq,2) ;
    end
end

nz = size(PI,1) ;
ns = size(PI,2) ;
nn = size(PHIn,2) ;

if numel(state)==nz
    state(:) = PI\state(:) ;
end

% check options / init
if ~isempty(state)
    assert(numel(state)==ns || numel(state)==nz,'Initial state vector must be consistent with state space system.') ;
else
    state = zeros(ns,1) ;
end
if isempty(shock)
    shock = zeros(nn,1) ;
end
if size(shock,2)==nn
    shock = shock' ;
elseif size(shock,1)==nn
    % excellent
else
    error('Shock vector/array must be conformable with state space system') ;
end
stp = abs(sum(shock,1))>eps ; % periods in which shocks occur
mps = [stp zeros(1,T-numel(stp))] ; % pad with zeros

dz = zeros(nz,T) ;
d2z = zeros(nz,T) ;
d3z = zeros(nz,T) ;

%% first order
q1 = zeros(ns,T) ;
% initial state
if ~isempty(state)
    q1(:,1) = state(:) ;
end
% response
for tt = 2:T
    if mps(tt-1)
        q1(:,tt) = PHIs*q1(:,tt-1) + PHIn*shock(:,tt-1) ;
    else
        q1(:,tt) = PHIs*q1(:,tt-1) ;
    end
    dz(:,tt) = PI*q1(:,tt) ;
end

if order>1
    %% second order
    q2 = zeros(nq2,T) ;
    % initial state
    if ~isempty(state)
        q2(:,1) = [zeros(ns,1); vech(state(:)*state(:)')] ;
    end
    % response
    for tt = 2:T
        if mps(tt-1)
            v2 = [vech(shock(:,tt-1)*shock(:,tt-1)'); vec(q1(:,1)*shock(:,tt-1)')] ;
            q2(:,tt) = T2q + T2qq*q2(:,tt-1) + T2qv*v2 ;
        else
            q2(:,tt) = T2q + T2qq*q2(:,tt-1) ;
        end
        d2z(:,tt) = T2z + T2zq*q2(:,tt) ;
    end
    
    if order>2
        %% third order
        q3 = zeros(nq3,T) ;
        % initial state
        if ~isempty(state)
            q3(:,1) = [zeros(ns,1); state(:); vec(state(:)*q2(:,1)')] ;
        end
        % response
        for tt = 2:T
            if mps(tt-1) 
                q3(:,tt) = T3q + T3qq*q3(:,tt-1) + T3qv*[shock(:,tt-1); vec(shock(:,tt-1)*q2(:,1)'); vec(q1(:,1)*v2'); vec(shock(:,tt-1)*v2')] ;
            else
                q3(:,tt) = T3q + T3qq*q3(:,tt-1) ;
            end
            d3z(:,tt) = T3z + T3zq*q3(:,tt) ;
        end
    end
    
    %% output
end

z = dz + (1/2)*d2z + (1/6)*d3z ;

end

